Local Self— Energy Approach For Electronic Structure Calculations 



N.E. Zein\ S. Y .Savrasov^, G. Kotliar^'* 
^RRC "Kurchatov Institute", Moscow, 123182, Russia 
^Department of Physics, University of California Davis, CA 95616, USA 
^Center for Material Theory, Department of Physics and Astronomy, 
Rutgers University, Piscataway, New Jersey, 08854, USA and 
'^Ecole Poly technique, 91128 Palaiseau Cedex, France 
(Dated: February 2, 2008) 

Using a novel self — consistent implementation of Hedin's GW perturbation theory we calculate 
space and energy dependent self-energy for a number of materials. We find it to be local in real 
space and rapidly convergent on second- to third- nearest neighbors. Corrections beyond GW are 
evaluated and shown to be completely localized within a single unit cell. This can be viewed as 
a fully self consistent implementation of the dynamical mean field theory for electronic structure 
calculations of real solids using a perturbative impurity solver. 
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The construction of a controlled practical approxima- 
tion to the many body problem of solid state physics is 
a long sought goal. Controlled approximations are im- 
portant because the accuracy of the results can be im- 
proved in a systematic way. This goal has been achieved 
in quantum chemistry by the configuration interaction 
(CI) method. CI can be thought as a controlled approx- 
imation that becomes more accurate as two factors are 
increased: a) the number of configurations (i.e. Slater 
determinants) kept and b) the size of the basis used to 
represent the one- particle orbitals which are used to 
represent the configurations. Dynamical mean field the- 
ory (DMFT) and its cluster extensions (C-DMFT) 
merge CI ideas with band structure methods. They allow 
us to tackle the problem of periodic infinite systems. 

The central goal is the computation of the one-electron 
Greens function, (its Fourier transform can be mea- 
sured via photoemission and inverse photoemission spec- 
troscopy), G{v,v' ,uj), and the self-energy I](r, r',a;). At 
the same time, following Hedin j^], one introduces the 
effective or screened interaction W{t,t' ,bj). The solu- 
tion of the full many-body problem can be formulated 
as the extremization of a functional L[G, W] = Tr In G + 
Ti\iiW~Y.G~IiW+<^[G, W]. It is defined as the Legen- 
dre transform of thermodynamic potential with respect 
to non-interacting Green function and bare Coulomb in- 
teraction N|. It strongly resembles the Luttinger-Ward 
functional 5] but has extremum both in the self-energy 
E and polarizability 11, which plays the role of the self- 
energy for W. 

The interaction functional $[G, W] is then expanded 
in a perturbative series. The first few graphs are shown 
in Fig^Ja) and corresponds to the Hartree, the GW, and 
the first correction beyond GW. Variations of $ over G 
and W give us S and 11. For the self-energy these dia- 
grams are given in Fig^^b). To solve the corresponding 
Dyson equations numerically one introduces a basis set 
and corresponding expansions for the self energies polar- 




FIG. 1: Diagrammatic representation of the interaction en- 
ergy functional $ (a) leading to the so called Hedin perturba- 
tion series for the self-energy (b). First term in the expansion 
is the Hartree contribution, second represents GW diagram, 
third is the correction to GW. 



izations and effective interactions. Cluster DMFT ideas 
truncate the functional <&, E and 11 by setting its vari- 
ables, i.e. the Greens functions, equal to zero beyond a 
given range R. When R is one lattice spacing we have 
the highly successful single site DMFT, as the range R 
increases the approximation converges to the solution of 
the full many-body problem. In this paper we address 
the central problem of determining the minimal range 
that is needed to obtain accurate results for various ma- 
terials. 

There are three different parameters that need to be 
increased to achieve convergence, a) the size of the basis 
set Lmax, b) the order of the perturbation theory kept 
n-max , c) the range of the graphs Rmax which needs to be 
kept to obtain accurate approximation. Rmax depends on 
L and n. We do not consider in this paper the important 
issue of convergence as a function of L„iax as well as the 
dependence of the range of the type of basis set chosen. 
Instead we make the choice of a minimal basis set and 
focus on the issues of convergence as a function of n and 
R. 
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Keeping Rmax equal to one lattice spacing and rimax = 
oo results in the Lmax orbitally degenerate single site 
DMFT approximation. Keeping Rmax — oo and n = 1 
corresponds to the famous GW approximation 0, 0, 0, 
H, I3- Rmax — 1 and Umax = 1 is rcduccd to the lo- 
cal GW approximation introduced by Zein and Antropov 
|lO|. as an approximation to accelerate the convergence 
of the GW method. Keeping Rmax{n = 1) = oo , 
Rmaxin > 1) = 1 and rimnsr. = oo constitutes the 
GW+DMFT approximationllfllilElIiai. Until now, 
this approximation has only been fully implemented in 
the context of model Hamiltonians [13 and its more re- 
alistic implementations jl2l | still contain adjustable pa- 
rameters such as the double counting correction. 

In this work we present a real-space cluster implemen- 
tation of the GW technique which allows us to monitor 
directly the locality of the self-energy in the real space. 
In addition, we evaluate corrections beyond the GW and 
answer the question of their convergency with respect to 
the cluster size. We show that the most non-local are the 
contributions from the diagrams with one loop. It was 
recognized early on^^l that the higher the order of the 
diagram, the more local it is because "crossing" integra- 
tion over internal wavevectors increases the role of large 
momentum leading to the locality in real space. We show 
here that in many real solids, the truncation of diagrams 
beyond one loop to the range of one lattice constant is 
already very accurate. 

We also investigate the smallest value of Rmax which 
is needed to obtain accurate results for each value of L. 
This allows us to obtain fully self consistent results, inde- 
pendent of the starting point where local density approx- 
imation (LDA) to density functional theory (DFT) 
serves in many cases |^ |^ . We establish that even 
in the case of semiconductors when the Coulomb inter- 
actions have an infinite range due to lack of screening, a 
reasonable small cluster produces very accurate results. 
The size of the cluster needed to obtain accurate results 
is reduced as the order in perturbation theory increases. 

For weakly correlated systems our approach can be re- 
garded as a trick to simplify and accelerate the solution 
of the GW equations and further perturbation correc- 
tions to it. Alternatively, our approach should be viewed 
as the first fully self-consistent implementation of an ab 



initio cluster DMFT method for solids (with second order 
perturbation theory playing the role of impurity solver). 

We discuss the results of our calculations for a num- 
ber of materials such as simple and transition metals as 
well as semiconductors. Especially the latter class rep- 
resents a hard case scenario for methodologies based on 
local self-energy approximations due to the long range 
nature of its statically screened Coulomb interaction. We 
focus on the electronic structure of Si, a benchmark in 
the past GW studies. Some aspects of this problem are 
still debated, such as the effect of higher lying energy 
states, core exchange, pseudopotential vs. all-electron 
approximations . We focus on the convergence of the 
electronic structure as a function of the cluster size used 
for determination of the self-energy for a given basis set. 
We also evaluate several diagrams beyond GW to exam- 
ine the issue of convergency of the whole perturbation 
theory with respect to the Coulomb interaction. 

Our implementation is based on the linear muffin-tin 
orbital (LMTO) method for electronic structure calcula- 
tions using the atomic sphere approximation (ASA) llTl 
which has been already used in realizations of GW [filial- 
The LMTO basis functions separate the wavevector and 
radial dependences 

X^(r) = $f(r)+^$^(r)5L- (1) 

L 

Here S^^ are the LMTO structure constants while ^^'"^ 
are linear combinations of the solutions of the radial 
Schrodinger equation as well as their energy derivatives 
taken at some fixed set of energies e^, at the center of in- 
terest which are matched continuously and differentiably 
to spherical Hankel (H) and Bessel (J) functions at the 
muffin-tin sphere boundaries. 

In order to find the matrix elements I]ct/3(k, icjn) of 
the self-energy operator E(r,r',iw„) in brackets of the 
LMTOs for a set of imaginary Matsubara frequencies 
ujn = {2n + 1)'kT it is useful to represent real space vec- 
tors r = p + H, r' = p' + R', where p, p' are restricted by 
the unit-cell, R, R' arc the lattice translations, and re- 
denote I](r, r', iujn) — Sfl(p, p', zcj„) (due to translational 
invariance we can always set R' = 0). Then, SQ,^(k, ia;„) 
has the following structure 



I]„^(k, za.„) = Si^^) (k, zc.„) + ^ J:if^ (k, ^^^)Slf, + ^ ^^^S^;^/' (k, zc.„) + ^ 5^ ^l;!^' (k, ^u;n)S^^ (2) 

L L L 

where E^'^?(k, ia;„) — J2 r ^''^^'^l^l^ ''■^n) and the cluster self-energy is given by the matrix element 

4't'(R,»tu„) = J^i^^*{p)J:B{p,p',^cu„)<^>[^){p')dpdp'. (3) 

I 

As we see, even if E(r, r', zc^n) can be local (i.e. non- zero only when both r and r' are in the same cell), the 
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TABLE L Correlational contribution to the self-energy ma- 
trix element, T,^Jf (R, iujo), in eV as a function of R for Fe, Ni, 
Na, Al and Si. Ri = a/y/2, R2 = a, R3 = ay/2 in fee struc- 
ture), Ri = ax/3/2 R2 = a, R3 = aV3/2 in bcc structure. 



Ri 


= 0.5a, R2 


= 0.83a R3 


= a in diamond structure. 








R = 


Ri 


R2 


R3 


Fe: 


S..(R) 


0.70 


0.35 


0.16 


0.00 


Fe: 


Edd(R) 


6.53 


0.05 


0.08 


0.00 


Ni: 


S.,(R) 


-0.54 


0.05 


0.03 


0.00 


Ni: 


Sdd(R) 


7.34 


0.38 


-0.03 


0.00 


Na: 


E,.(R) 


-1.36 


0.52 


0.30 


0.05 


Al: 


Epp(R) 


0.46 


0.16 


-0.08 


0.00 


Si: 


S,.(R) 


0.95 


0.30 


0.14 


0.03 


Si: 


Epp(R) 


-1.06 


0.05 


-0.27 


-0.03 



matrix elements T,ap(k.,iujn), Eq. acquire some k- 
dependence through the structure constants. It is due to 
tails of the basis functions extended over all space, and it 
is quite analogous to the k-dependence of local poten- 
tial matrix elements in LDA. Such k-dependence can be 
called " kinematical" , and in the following we will distin- 
guish it from the dynamical k dependence connected to 
the existence of R ^ elements of p', zwn) which 

is the main focus of the present work. 

From Eq. Q follows that the one-electron Green func- 
tion can be represented in a factorized form 

fj,iy—HJ 

(4) 

and, as a result, the polarization operator n(r, r',iw„), 

has a similar structure. We solve the equation for the 
dynamically screened interaction W — V — VHW on 
the product of basis functions following Ref. [H Af- 
ter finding M^(r,r',iw„) the self-energy is calculated ei- 
ther as S(R) — G(R)VF(R) (R-space version) or as 
= J2kG{k)W{k + q) (k-space version). All calcu- 
lations are performed on the imaginary axis. Due to large 
frequency behavior of the Green function proportional to 
l/icOn a special care should be taken of the direct ex- 
change contribution to the self-energy (E^, = —GV ), as 
the sum over large u;„ needs to be done analytically. The 
remaining portion, Sf, = G{VIIW), is due to correlations, 
and the sum over internal frequencies is rapidly conver- 
gent. Finally, in order to obtain the electronic spectrum 
for real frequencies we analytically continue the Green 
function using Pade's approximation. 

We perform self-consistent GW calculations for Fe, Ni, 
Na, and Si using our newly developed cluster algorithm 
and obtain the self-energies in real space. To discuss the 
results of these calculations, TablelUlists the diagonal ma- 
trix elements of the correlational part of the self-energy 
^cTi^'' ("l^' *^o), (0, as a function of R for the value of 



TABLE II: Polarizability nij(R,aj = 0) as a function of clus- 
ter size R as well as Ilij (k — 0,ij = 0) for Si {ij— 1,2 numerates 
Si atoms, jj=3,4 numerates empty spheres) . Yl n shows how 
the sum rule n(A; 0, a; = 0) ^ is fulfilled in both R-space 
and k-space calculations. 







ni2 


ni3 


ni4 


En 


R = 


5.34 











14.78 


Ri 


4.98 


-2.60 


-0.55 


0.21 


6.76 


R2 


4.98 


-2.67 


-1.00 


-0.09 


3.45 


R3 


4.98 


-3.33 


-1.159 


-0.21 


0.81 


R4 


4.96 


-3.36 


-.14 


-0.23 


0.41 


k-space 


4.88 


-3.40 


-1.21 


-0.26 


10-" 



ujQ = ttT = tt400K. (In case of magnetic ground state 
Y.^^{iujo) is for majority spins.) From Table ^ it fol- 
lows that Ec,rf(i(R = 0) dominates in transition metals 
and falls off very quickly in nearest neighbors, because of 
both small overlap between d-orbitals and large screen- 
ing at small energies. The Sc,ss(R-) in these metals falls 
off more gradually but its value is negligible in compari- 
son with d. It explains why recently developed one-site 
approach ^3 is successful in this case. In simple metals 
like Na and Al only and Sp are significant, they fall 
off gradually, but nevertheless on the third coordination 
sphere they become small as compared to E(R = 0), even 
if they are multiplied by a number of nearest neighbors 
at this sphere. In Si, both Ec,ss(R) and Sc,pp(R) fall off 
pretty slowly but nevertheless become very small at R's 
comparable with the size of the unit cell in accord with 
the conclusions reached in the pioneering work jS| • 

A proper approximation of the polarizability 11 of a 
semiconductor, requires going beyond the local approx- 
imation in order to fulfill the f sum rule n(k — > 0, w = 
0) 0. This requires cancellations between local and 
non-local terms. In Tabled we show how this sum rule 
is fulfilled for Si, where values of ny(R, w = 0) {ij nu- 
merate atoms in the unit cell) are listed for several R's. 
As we utilize the atomic sphere approximation, in order 
to reach close packing we consider two Si atoms (ij=l,2) 
and two empty spheres (ij=3,4). Table ITU also lists the 
matrix elements ny(k = 0,a; = 0) obtained by the k- 
space calculation as well as sums, ^ Ft, over the atomic 
coordinates, which show how the sum rule is fulfilled. We 
see that while the k -space version automatically leads to 
obeying the sum rule, the R -space version needs cluster 
sizes extended up to four nearest neighbors. 

We now investigate the convergency of the electronic 
structure of Si using the cluster GW method. Both the 
direct and indirect energy gaps as well as the valence 
band width are calculated by varying the size of the clus- 
ter used when evaluating correlational part of the self- 
energy. The behavior of the electronic structure as a 
function of the cluster size is schematically shown on Fig 
121 They indicate that the non-locality of the self-energy 
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TABLE III: Comparison between calculated energy gaps (eV) 
in semiconductors and insualtors using non-self-consistent 
"first shot" (FS) and self-consistent (SC) GW methods ob- 
tained in the present work, results of other available GW cal- 
culations and experiment. 



0.5 0.83 1 R=oo 
Cutoff Radius R 



FIG. 2: Direct and indirect energy gaps as well as valence 
bandwidth for Si in eV calculated by varying the size of the 
cluster used when evaluating correlational part of the self- 
energy within the GW method. 



in Si approximately extends up to a third coordination 
sphere. 

We have evaluated the first self-energy correction be- 
yond GW, AI](R), for Si as a function of the cluster size. 
Our obtained on-site values for s and p electrons are as 
follows: AEs(O) =0.22 eV and AEp(O) =-0.14 eV. The 
self-energies at first and second coordination sphere are 
given by: ASs(Ri) =0.002 eV, AEp(Ri) =-0.002 eV, 
AEs(R2) =0.0003 eV, AEp(R2) =-0.0005 eV. We see 
that the correction to GW is completely local and per- 
mits us to calculate it in the real space vifhich much less 
time-consuming. It is also interesting to note that the 
energy dependence of this correction exists but at the 
same scale as in 'Saw- 

We finally would like to address a highly interesting 
question on the accuracy of the GW calculations in pre- 
dicting the energy gaps in various semiconductors and 
insulators. We have performed such calculations with 
and without imposing the self-consistency for the Green 
function which in the latter case corresponds to the LDA 
Green function used in evaluating the interactions and 
the self-energies. We refer to these calculations as self- 
consistent (SC) and "first shot" (FS) ones. We include 
the core exchange effects, whose importance has been 
recently pointed out 0. The results of these studies 
are presented in Table mil where we list the obtained en- 
ergy gaps for a whole series of materials such as C, Si, 
MgO, and AlAs. The available previous calculations and 
experimental data are also listed in this table for com- 
parison. Despite the possible inaccuracies introduced by 
the use of a limited basis set, which does not include very 
high-energy states, we see that our non-self-consistent 
calculations are in good agreement with the published 
data. The self-consistency is not only of conceptual im- 
portance, since it makes the solution of the problem in- 





FS (this work) 


SC (this work) 


FS 


SC 


Exp. 


c 


5.00 


5.02 


4.92" 




5.48** 


Si 


0.86 


1.10 


0.85" 


1.05" 


1.17' 


MgO 


8.00 


5.90 


8.3'' 




7.8'' 


AlAs 


1.33 


1.90 


1.05" 




2.24' 



"Rcfl;'Rcf 20;'=Rcfi6; ''Ref. IH'^Reflll; 



dependent on the starting point, but can also affect the 
quality of the results in some systems such as AlAs. 

In conclusion we have developed a self-consistent clus- 
ter DMFT methodology which allows us to monitor the 
locality of the self-energy in the real space. As a first 
application of the method we evaluated first contribu- 
tion in Hedin's perturbation series for the self-energy in 
Si and found it to be completely local. Our calculated 
energy gap values for a number of semiconductors and 
insulators are found in good agreement with experiment. 
The fast convergence in real space demonstrated in this 
paper simplifies and accelerates the GW approximation 
which is useful in many solids. Our approach is also the 
ideal starting point for resuming the perturbation theory 
at the local level using non-perturbative solvers as was 
done in Ref. This step will allow us to study very 

strongly correlated systems completing the first princi- 
ples C-DMFT program for solids. 
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